Coexistence of regular and irregular dynamics in complex networks of pulse-coupled 

oscillators 
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For general networks of pulse-coupled oscillators, including regular, random, and more complex 
networks, we develop an exact stability analysis of synchronous states. As opposed to conventional 
stability analysis, here stability is determined by a multitude of linear operators. We treat this multi- 
operator problem exactly and show that for inhibitory interactions the synchronous state is stable, 
independent of the parameters and the network connectivity. In randomly connected networks with 
strong interactions this synchronous state, displaying regular dynamics, coexists with a balanced 
state that exhibits irregular dynamics such that external signals may switch the network between 
qualitatively distinct states. 

PACS numbers: 89.75.-k, 05.45.-a, 87.10. +e 
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Complex networks appear as a variety of natural and 
artificial systems, ranging from the world wide web and 
electrical power grids to metabolic and neural networks 
[HQ- While recent studies have focussed on their struc- 
ture pj, the dynamics in such networks constitute a chal- 
lenging issue of current and future research Even 
if the individual vertices of the network are simple dy- 
namical systems, such as limit cycle oscillators, an ex- 
act mathematical analysis of their collective dynamics is 
often highly intricate, due to the complex connectivity 
structure. 

As a prototypical class of dynamical systems interact- 
ing on networks, pulse-coupled units have received a sig- 
nificant amount of interest because of their relevance to 
diverse natural systems 0. QJ. 0- 0- 01 including cardiac 
pacemaker cells, flashing fireflies, earthquakes, and bi- 
ological neural networks. Particularly in neuroscience, 
these models |g] are essential for understanding collec- 
tive dynamic phenomena such as synchronization or the 
propagation of sensory signals through extended net- 
works |S|,Ui|. Although biological neural networks, like 
other networks occurring in the real world, often possess a 
complex connectivity structure, most theoretical studies 
on pulse-coupled oscillators are either restricted to net- 
works of globally coupled oscillators and simple regular 
networks, or work in some mean field limit 0- Q- 0- 10| - 

In this Letter, we study pulse-coupled oscillators in- 
teracting on networks with general connectivities, includ- 
ing fully connected, regular, random, and more complex 
topologies. We develop an exact stability analysis of the 
synchronous state. In contrast to conventional stabil- 
ity problems, the first order stability operator here is 
not linear, but can be expressed by a multitude of dis- 
tinct linear operators, the domains of which depend on 
the rank order of a specific perturbation. For generally 
structured networks, the number of operators increases 
exponentially with the size of the network. Our analysis 
provides a method to treat this multi-operator problem 
analytically. For networks with inhibitory couplings, we 
prove that the synchronous state is stable, independent 



of the parameters and the connectivity structure. 

Proceeding from this result, we show that the syn- 
chronous state, that displays regular dynamics, coexists 
with a state of highly irregular dynamics in randomly 
connected networks. We suggest a simple mechanism for 
switching between these states. These results establish 
that the behavior of networks of pulse-coupled units for 
a given set of parameters may be dominated by qualita- 
tively distinct dynamical states. 

We consider a system of N coupled oscillators Q which 
interact on directed graphs by sending and receiving 
pulses. The structure of this graph is specified by the 
sets Pre(i) of presynaptic oscillators that send pulses to 
oscillator i. These sets determine the sets Post(i) of post- 
synaptic oscillators that receive pulses from i. A phase- 
like variable <j>i{t) G (— oo, 1] specifies the state of each 
oscillator at time t. The dynamics of a non-interacting 
oscillator i is given by 



i(t)/dt = 1. 



(1) 



When oscillator i reaches a threshold, <f>i(t) = 1, its phase 
is reset to zero, <j>i{t + ) = 0, and the oscillator is said 
to "fire". A pulse is sent to all postsynaptic oscillators 
j E Post(i) which receive this signal after a delay time 
r. Depending on whether the input is subthreshold or 
suprathreshold, the incoming signal induces a phase jump 



&((t + T)+) :=min{C/- 1 (C/(^(t + r))+£, l ),l} (2) 

which depends on the instantaneous phase <pj (t + r) of 
the postsynaptic oscillator and the coupling strength ej%. 
The phase dependence is determined by a twice con- 
tinuously differentiable 'potential' function U((j>) that is 
assumed to be strictly increasing, U'{<f>) > 0, concave 
(down), U" ((/>)< 0, and normalized such that £7(0) = 0, 
U(l) = 1 (cf. i). 

By choosing an appropriate function U, this model is 
equivalent (cf. |.]jj) to different well known models of in- 
teracting threshold elements. For instance, for the leaky 



2 



integrate-and-fire oscillator defined by the linear differen- 
tial equation V = I — V ^| (and threshold at V = 1), one 
obtains U IF ((f>) = I(l-e LT ^) where Tip = log(7/(I-l)) 
is the period of a non-interacting oscillator and I > 1 is 
a suprathreshold external current. Oscillators described 
by nonlinear differential equations are covered by the 
Mirollo-Strogatz approach, too. For instance, the con- 
ductance based threshold model of a neuron |3] leads to 
a different, more complicated function Ucb(4) (for de- 
tails see All analytical results presented here are 
derived for the above, general class of interaction func- 
tions. In numerical investigations, we use the functional 
form Uif but find qualitatively similar results for differ- 
ent U. In this Letter, we focus on inhibitorily coupled 
networks (all < and en = 0). 

We perform a stability analysis of the synchronous 
state (4>i(t) — 4>o{t) for all i) that exists if the coupling 
strengths are normalized such that X^ePrc(i) £ ij = £ < 0. 
Its period is given by 

T = t + 1 - /? (3) 

where 0o = U~ 1 (U(t) + e). To construct a strobo- 
scope map, a perturbation 8(0) = 8 = (8%, . . . ,8n) of 
the phases, defined by 



&(0) = mo) + s 1 , 



(4) 



is ordered according to the rank order rank(<5) of the 8n 
For each oscillator i we label the perturbations 5j of its 
presynaptic oscillators j (for which £y ^ 0) according to 
their size 



Ai,i > A h2 >■■■> A lM 
Pre(i)| is the number of its 



(5) 



)resynaptic 
1211 . In ad- 



where ki :— 

oscillators, called in-degree in graph theory 
dition, we define A^o = 8$. For illustration assume that 
an oscillator i has exactly two presynaptic oscillators ji 
and j2 such that Pre(i) = {31,32} and ki — 2. For cer- 
tain perturbations, oscillator i first receives a signal from 
oscillator 32 and slightly later from oscillator ji . This de- 
termines the rank order (8j 2 > 8j ± ) such that A^i = Sj 2 
and Aj t 2 = 8j t . 

Using the phase shift function h ((p, e) := U^ 1 (U((j))+e) 
and denoting Di iTl := Aj. n _i — Aj >n for n G {1, . . . , ki} we 
compute the time evolution of phase-perturbations Si <C 
1, starting near </>o(0) = r/2 without loss of generality. 
The stroboscopic time-T map of the perturbations, Si 1— » 
Si(T), is obtained from the scheme 
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where the right column gives the phases <fo (t) of oscillator 
i at times t of pulse receptions or reset given in the left 
column. Here the presynaptic oscillator from which oscil- 
lator i receives the n th pulse during this cycle is labeled 
by j n . The time to threshold 
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- - A i>fei + 1 - (3 tM 



(6) 



is always smaller than the period T. Hence the period- T 
map of the perturbation can be expressed as 

Si (T) = T- 7f } - I = f3 iM - ft, + A itki . (7) 

Expanding 0^ for small Di^ n one can prove by induction 
that to first order 



0i,ki = 00 + y^Pi,n-lA.-, 



(8) 



where 



P " l -~ U'(U-i(U(t) + e)) [) 
for n G {0, 1, . . . , ki}. This results in a first order map 

8(T) = AS (10) 
where the elements of the matrix A are given by 

G Pre(i) 



Pi,n ~ Pi, n -1 if J = jr. 
An = < Pi,o if i = i 

if 3 i Pre(i) U {i}. 



(11) 



Since j n in © and (1 1 H identifies the n th pulse received 
during this cycle by oscillator i, the first order oper- 
ator depends on the rank order of the perturbations, 
A = A(rank(8)), and the map A8 is piecewise linear. 
In general, signals arriving almost simultaneously at the 
same oscillator induce different phase changes, depending 
on the order of arrival: For the above example of an os- 
cillator i with exactly two presynaptic oscillators ji and 
32 and equal coupling strengths, Eij 1 = Sij 2 , the first 
of the two arriving signals has a larger effect, encoded 
in the pi iU , than the second, by virtue of the concavity 
of U(4>) (cf. Eq. (2J). The respective matrix elements 
Aij 1 and Aij 2 are differences between certain Pi >n and 
therefore have different values depending on which sig- 
nal is received first. This is induced by the structure 
of the network together with the jump-like interactions. 
For networks with homogeneous, global coupling differ- 
ent matrices A can be identified by appropriately per- 
muting the oscillator indices. In general, however, this is 
impossible. 

Hence, in this stability problem, given a network struc- 
ture, one generally has to deal with an exponential num- 
ber of operators instead of a single stability matrix. We 



treat all these operators simultaneously: It is straight- 
forward to show that for all matrices A (independent of 
the rank order of a perturbation and the parameters) 
the matrix elements are non-negative, Ay > 0. Due 
to time-translation invariance all A are normalized row- 
wise, J2j Aij — 1 for all i, and exhibit a trivial eigenvalue 
Ai = 1. Moreover, the diagonal elements are identical 
and smaller than one, An = Ao < 1. The synchronous 
state is thus stable, because the inequality 

maxi | Si(T) | < max* Y^j I Aj II^j I 

< maxi Y,j \ A %j I max fc 141 = max fc \5 k \, 

(12) 

is satisfied for all matrices A. 

For a convex potential function U, where U" > 0, and 
excitatory interactions (ey > 0) the synchronous state is 
stable as well. The above proof applies if the total input 
e is not suprathreshold, i.e. e < 1 — U(t). Thus whereas 
excitatory interactions must not be too strong for ap- 
plicability of the proof, inhibititory interactions may be 
arbitrarily strong. 

For structural stability of the stable synchronous state 
it is required that the non-trivial eigenvalues of the ma- 
trices A are separated from the unit circle. An instruc- 
tive example is given by a network of integrate-and-fire 
oscillators, U{<j>) — Uxf(4>), where all matrices l(TT|l are 
degenerate if £y = e/fcj for all j S Pre(i). In this 
case, the eigenvalues of a single matrix completely char- 
acterize the dynamics in the vicinity of the synchronous 
state. Numerically, we find that in large random net- 
works in which every connection is present with proba- 
bility p all non-trivial eigenvalues are located in a disk 
D = {z G C| \z — Aq\ < r} of radius r that is centered at 
Aq < 1 and separated from the unit circle. An estimate 
for the radius, r = (1 - A )(p" 1 - 1) 1/2 ^ 1/2 for N > 1, 
can be obtained from the theory of Gaussian asymmet- 
ric random matrices 0] . We find that this estimate well 
agrees with our numerical results [3. This indicates 
that in the limit of large N, all non-trivial eigenvalues 
are concentrated near z = Aq and thus separated from 
unit circle. 

The above analysis shows that for inhibitory coupling 
the synchronous state is stable, independent of the pa- 
rameters and the network structure. Numerical simula- 
tions show that for a network at given parameters this 
synchronous state often coexists with one or more other 
attractors. A particularly important example which oc- 
curs in randomly connected networks with strong inter- 
actions, is a balanced state (cf. 0,^E|) that exhibits irreg- 
ular dynamics. In this balanced state, found originally in 
binary neural networks > inhibitory and excitatory in- 
puts cancel each other on average but fluctuations lead to 
a variability of the membrane potential and a high irreg- 
ularity in firing times (see also Q). Figures QJi,b display 
sample trajectories of the potentials U(<f>i) of three oscil- 
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PIG. 1: Coexistence of (a) synchronous and (b) irregular dy- 
namics in a random network (N = 400, p = 0.2, / = 4.0, 
e = 16.0, r = 0.035). (a),(b): Trajectories of the poten- 
tial U{4>i) of three oscillators (angular bars: time scale (hor- 
izontal) At = 8; potential scales (vertical) (a) AU = 8, (b) 
AU = 2 ; spikes of height AU = 1 added at firing times). 
(c),(d): Distributions (c) p v of rates and (d) pcv of the co- 
efficient of variation, displayed for the irregular (dark gray) 
and synchronous (light gray) dynamics. 

lators for the same random network, making obvious the 
two distinct kinds of coexisting dynamics. 

The dynamical differences are quantified by a his- 
togram p v of oscillator rates (Fig. Qt) 

v i = ((ti,n+l — ti,n) n ) , (13) 

the reciprocal values of the time averaged inter-spike- 
intervals. Here the ti >n are the times when oscillator 
i fires the n th time. The temporal irregularity of the 
firing-sequence of single oscillators i is measured by the 
coefficient of variation 

CV, = {v} ({U. n+1 - U. n f) n - I)* , (14) 

defined as the ratio of the standard deviation of the inter- 
spike intervals and their average. A histogram pcv of 
the CVi (Fig. QJl) shows that the irregular state exhibits 
coefficients of variation near one, the coefficient of vari- 
ation of a Poisson process. Such irregular states occur 
robustly when changing parameters and network topol- 
ogy; on the other hand, the size of the basin of attraction 
of the synchronous state is also significant and increases 
with increasing delay r . 

The coexistence of two qualitatively different kinds of 
dynamics leads to the question how regular dynamics 
can be induced when the system currently is in an ir- 
regular state and vice versa. A simple mechanism to 
synchronize oscillators that are in a state of irregular 
firing is the delivery of two sufficiently strong external 
excitatory (phase-advancing) pulses that are separated 
by a time At <G (r, 1), cf. Fig. El The first pulse then 
leads to a synchronization of phases due to simultaneous 
suprathreshold input (cf. Eq. (01) . If there are traveling 
signals that have been sent but not received at the time 
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addressing the dynamics in networks, now that impor- 
tant aspects of their complex structure have been under- 
stood m 
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FIG. 2: Switching between synchronous and irregular dynam- 
ics (N = 400, p = 0.2, I = 4.0, e = 16.0, r = 0.14). Firing 
times of five oscillators are shown in a time window At = 240. 
Vertical dashed lines mark external perturbations: (i) large 
excitatory pulses lead to synchronous state, (ii) a small ran- 
dom perturbation (|A</>i| < 0.18) is restored (iii) a sufficiently 
large random perturbation (|A</>j| < 0.36) leads to an irregu- 
lar state. Bottom: Time evolution of the spread of the spike 
times after perturbation (ii), total length At — 0.25 each. 



of the first pulse, a second pulse after a time At > t 
is needed that synchronizes the phases after all inter- 
nal signals have been received. This synchronous state 
is not affected by small random perturbations, whereas 
large random perturbations lead back to irregular dynam- 
ics (Fig. EJ. Mechanisms for both directions of switching 
may be realized in biological neural networks by exter- 
nal neuronal populations: While strong external pulses 
may be generated by external neurons that are highly 
synchronized, a random perturbation can be realized by 
neurons which fire irregularly. 

Most previous studies of the dynamics of networks 
of pulse-coupled units focussed on regular networks or 
worked in some mean field limit These studies 

often relied on the analysis of bifurcations from one state 
to another as an external parameter is changed. Based on 
the stability analysis developed here, that applies to net- 
works with general connectivity, we have demonstrated 
that regular synchronous dynamics may coexist with ir- 
regular dynamics in sufficiently complex networks. The 
coexistence of qualitatively different states at identical 
parameters indicates that bifurcation approaches may of- 
ten not give a complete picture of the network dynamics, 
if the network structure is too complex. This fact may 
well apply not only to networks of pulse-coupled units but 
also to the dynamics of many other complex networks. In 
addition, our results emphasize that in complex networks 
of pulse-coupled units the occurrence of temporally regu- 
lar and irregular firing patterns may typically reflect the 
collective state of the network rather than the dynamics 
of individual units. 

The analysis presented in this Letter demonstrates 
that the dynamics in certain complex networks can be 
revealed by considering the vertices as units with simple 
dynamical properties, e.g. intrinsic oscillators. Such sys- 
tems provide promising starting points for future studies 
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